Mean-field description of ground-state properties 

of drip-line nuclei. 
(I) Shell-correction method 



W. Nazarewicz, a ' 6,c,d , T.R. Werner, a,b,c and J. Dobaczewski, c 
a Physics Division, Oak Ridge National Laboratory, 
P.O. Box 2008, Oak Ridge, Tennessee 31831, USA 

b Department of Physics and Astronomy, University of Tennessee, 
Knoxville, Tennessee 37996, USA 

c Institute of Theoretical Physics, Warsaw University, 
ul. Hoza 69, 00-681 Warsaw, Poland, 

d Joint Institute for Heavy-Ion Research, Oak Ridge, Tennessee 31831, USA 



Abstract 
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I. INTRODUCTION 



The theoretical description of drip-line nuclei is one of the most exciting challenges today. 
The coupling between bound states and continuum invites strong interplay between various 
aspects of nuclear structure and reaction theory, and calculations have strong astrophysical 
implications, especially in the context of the r-process mechanism 

Theoretically, the physics of nuclei with very large values of neutron or proton excess 
is a challenge for well-established models of nuclear structure and, because of dramatic 
extrapolations involved, it invites a variety of theoretical approaches. Since the parameters 
of interactions used in the usual shell-model or mean-field calculations are determined so as 
to reproduce the properties of beta-stable nuclei, the parameters may not always be proper 
to be used in the calculation of drip-line nuclei. One hopes, however, that spectroscopy of 
exotic nuclei by means of a radioactive ion beam technique will lead to a better determination 
of forces, at least those interaction components that depend on isospin degrees of freedom. 

The closeness of the Fermi level to the particle continuum makes the theoretical de- 
scription of drip-line nuclei a very challenging task. To put things in some perspective, 
Fig. [I] displays the average potential wells, characteristic of a typical beta-stable system 
( 120 Sn), a neutron drip-line system ( 150 Sn), and a proton drip-line system ( 100 Sn). While 
the low-energy structure of 120 Sn is almost exclusively determined by the particle-hole or 
pair excitations across the Fermi level from bound states to bound states (corrected for po- 
larization effects due to giant vibrations), the lowest particle-hole or pair modes in drip-line 
nuclei are embedded in the particle continuum. Consequently, any tool of nuclear structure 
theory that aims at describing many-body correlations starting from the mean-field-based 
single-particle basis (such as shell-model, BCS, RPA, etc.) has to be modified in the new 
regime. 

One of the most important nuclear properties is its mass. The ability of a theoretical 
model to reproduce the nuclear binding energy is its ultimate test; it determines its reliability 
and practical usefulness. There exist many mass calculations (e.g., those based on the nuclear 
shell model) focused on a narrow region of the nuclear chart. These calculations are very 
successful in reproducing the data in a given region, but their applicability to other nuclei 
is limited. In order to extrapolate far from stability, the large-scale global mass calculations 
are usually used (see, e.g., reviews |2|U). Since their parameters are optimized to reproduce 
known atomic masses, it is by no means obvious whether the particle number dependence 
predicted by global calculations at very large (or very small) values of the relative neutron 
excess 



is correct. In this context, a good example is the inclusion of the Coulomb redistribution 
energy term in the finite-range droplet model - strongly motivated by the recent mass 
measurements for the heaviest elements [Q. 

As far as the placement of nuclear drip lines is concerned, it is not the absolute value 
of nuclear mass but rather the mass difference between two isotopes that is of interest. 
The difference between binding energies, B(Z, N), determines both the one- neutron separa- 
tion energy, S n (Z, N) = B n (Z, N) — B n (Z, N — 1), and the two-neutron separation energy, 
S2n{Z, N) = B n (Z, N) — B n (Z, N — 2). If S n becomes negative, then a nucleus is neutron 
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unstable; the condition S n =0 determines the position of the one-neutron drip line. By the 
same token, if becomes negative, a nucleus is unstable against the emission of two neu- 
trons; the condition S , 2„=0 determines the position of the two-neutron drip line. (Proton 
drip lines are defined analogously.) Since S n and S 2n are differences between two masses, 
the corresponding statistical rms deviations are greater than the estimated rms deviation for 
the total binding energy. That is, these quantities are less certain than calculated masses. 

Because of their sensitivity to various theoretical details, predicted drip lines are strongly 
model-dependent The variations between predictions can be attributed to: (i) fun- 

damental differences between microscopic approaches, (ii) different effective interactions 
employed within the same approach, and (iii) different approximations used when solving 
the nuclear many-body problem within a given approach. In this study, several aspects of 
nuclear structure at the limits of extreme isospin are discussed by means of the macroscopic- 
microscopic approach. To test the influence of the particle continuum on shell corrections 
and pairing energies, we use the two versions of the shell-correction method with the Woods- 
Saxon average potential (first version is based on the standard averaging method and the 
second version is based on the semiclassical Wigner-Kirkwood expansion). In Part II of 
this work | 10| , pairing properties of drip-line systems and the sensitivity of predictions to 
effective forces are investigated by means of self-consistent Hartree-Fock (HF) and Hartree- 
Fock-Bogolyubov (HFB) approaches. 



II. SHELL-CORRECTION METHOD AND PARTICLE CONTINUUM 



The main assumption of the shell-correction (macroscopic-microscopic) method [pTT| [15 
is that the total energy of a nucleus can be decomposed into two parts: 

E = E + -Eshell, (2.1) 

where E is the macroscopic energy (smoothly depending on the number of nucleons and 
thus associated with the "uniform" distribution of single-particle orbitals), and -E s heii is the 
shell-correction term fluctuating with particle number reflecting the non-uniformities of the 
single-particle level distribution, i.e., shell effects. In order to make a separation ( [2.1|) , one 
starts from the one-body HF density matrix p 

p(x', x)=Y, <f>i(x')<l>*(x), x = (r, <r) , (2.2) 

i—occ 

which can be decomposed into a "smoothed" density p and a correction 5p, which fluctuates 
with the shell filling 

p = p + 5p. (2.3) 
The density matrix p can be expressed by means of the smoothed distribution numbers n« 

p(x',x) = J2nMx')<j ) *(x). (2.4) 

i 

When considered as a function of the single-particle energies e,, the numbers hi vary smoothly 
in an energy interval of the order of the energy difference between major shells. The ?Vs 
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can, in principle, take values smaller than zero and larger than unity (T^j. Consequently, 
they do not represent the single-particle occupations in the strict sense. 

In order to justify Eq. ( [2.1Q , the expectation value of a HF Hamiltonian (containing the 
kinetic energy, t, and the two-body interaction, v) can be written in terms of p and Sp 
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E UF = Tr(tp) + ^Tr Ti(pvp) = E + E osc + 0(5p 2 ), (2.5) 



where 



E = Tr(tp) + ^Tr Tr(pvp) (2.6) 

is the average part of E'hf and 

E osc = Tr(h HF 6p) (huF = t + Tr{vp)) (2.7) 

is the first-order term in Sp representing the shell-correction contribution to £7hf- If a 
deformed average potential gives a similar spectrum to the averaged HF potential /ihf then 
the oscillatory part of -Ehf, given by Eq. (|2.7|) , is very close to that of the deformed shell 
model, E she n^iE osc . The second-order term in Eq. (|2.5| ) is usually very small and can be 



neglected [|17]]. The above relation, known as the Strutinsky energy theorem, makes it possible 
to calculate the total energy using the non-self-consistent, deformed independent-particle 
model; the average part E is usually replaced by the corresponding phenomenological liquid- 
drop (or droplet) model value, E mauCT . It is important that -E s heii must not contain any regular 
(smooth) terms analogous to those already included in the phenomenological macroscopic 
part. (The extension of the energy theorem to the case with pairing is straightforward. 
The resulting expression for shell correction contains an additional contribution from the 
fluctuating part of the pairing energy, see Sec. |V D| .) 

There are two single-particle level densities that define the shell correction. The (de- 
formed) shell-model single-particle level density 

g(e) = Y,5(e-e i ) (2.8) 

i 

gives the single-particle energy, E sp ,. The smooth single-particle energy, E sp ,, is given by 
the mean single-particle level density, g(e), obtained from g(e) by folding with a smoothing 
function f(x): 

m = - r deW)f (—) =-E/ (—) ■ ( 2 -9) 



\ 7 

In Eq. ( |2.9| ), 7 is the smoothing range; it should be larger than the typical distance between 
major shells. As follows from Eqs. (|2.7|) - (|2.9|) , the shell correction can be calculated by 
taking the difference between the sum of occupied levels and its average value, i.e. 

£sheii = E s . p . - E s . p . = *i - t ^(e)de, (2.10) 

i-occ J -°° 
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where A is the smoothed Fermi level defined through the particle number equation: 

N= [ X ~g(e)de. (2.11) 



The folding function f(x) can be written as a product 

f(x)=u(x)P p (x), (2.12) 

where uj(x) is a weighting function and P p (x) is the so-called curvature-correction polynomial 
of the pth order. The smoothing procedure should be unambiguous, i.e., the averaging should 
extract only the fluctuating part, leaving the smooth part untouched. This condition defines 
P p (x) for any specific choice of u(x) . For instance, for the infinite potential well and a 
Gaussian weighting function, lo(x) = 7r _1 ' 2 exp(— x 2 ), the curvature-correction polynomial 
can be be expanded in a finite series of Hermite polynomials of even order 0,0] 

p (-l) k / 2 

Pp{x) = k ^ 2 W/iy Hk{x) - (2,13) 

(For other examples of f(x), see Ref. [PJ.) 



The smoothed single-particle energy can be expressed in the form |L9] : 

M dE 
E s . p .= e ~g(e)de = Y,e l n l + l^, (2.14) 

J-oo Y "7 

where the smoothed distribution numbers are 

fk = - t */f— V ( 2 - 15 ) 



7 J-oo \ 7 

Since the value of E s . v , should not depend on the smoothing range 7 (nor on the order of 
curvature correction p), the second term in Eq. ( [2.14|) must vanish, i.e. 

^ = (^h = oV (2.16) 
07 \ dp J 

If the above plateau condition does not hold, the Strutinsky averaging method does not yield 
an unambiguous result. 

The treatment of the particle continuum in the shell-correction approach is an old prob- 
lem. For finite-depth potentials, the sum in Eq. ( |2.9| ) should be replaced by a sum over 
bound single-particle states and an integral over positive-energy single-particle continuum, 
i.e. 

E +J deA 9coilt (e), (2.17) 

i i— bound " 



where Ag cont is the continuum shell level density 120,21 



A<w(e) = (2.18) 
7r ae 
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and S e is the corresponding phase shift (calculated for all partial waves). 

Already in 1970 Lin [22j pointed out that the contribution from the particle continuum 



can affect the value of .E'sheib even for nuclei at the beta-stability line. In particular, he 
considered the positive-energy continuum in the energy interval < e < 20 MeV to calculate 
the neutron shell correction for 208 Pb. No plateau in E s . p . was obtained and the result 
turned out be strongly 7- and p-dependent. Soon afterwards Ross and Bhaduri |23| carried 



out calculations along the same lines as in Ref. \12\ and demonstrated that, by taking into 
account contributions from all neutron resonances up to ~100 MeV in 208 Pb, the plateau 
condition (|2.16|) could be met. 



Bolsterli et al. made an attempt to simulate the effect of the continuum by using 
the quasibound states, i.e., the states resulting from the diagonalization of a finite potential 
in a large harmonic oscillator basis. The authors suggested a "working" prescription (the 
number of oscillator shells should not be too large, iV osc ~12), but no systematic analysis was 
carried out. (Later it was shown in Ref. [^3J] that the inclusion of quasibound states requires 
relatively high correction orders.) In the following, the Strutinsky averaging procedure 
including quasibound states will be referred to as the standard averaging method (SAM). 

There have been some suggestions on how to generalize the Strutinsky averaging proce- 
dure for finite potentials. Bunatian et al. [|TJ]] exploited only the bound states of a finite-depth 
single-particle potential and derived modified expressions for the curvature correction. This 
method was then improved and exploited in series of papers by Strutinsky and Ivanyuk 
|25|-p7 . For instance, Ref. |27j demonstrates that the standard way of computing shell 



corrections, based on quasibound states, leads to serious errors in theoretical predictions 
for masses, fission barriers, and deformation energies. They also pointed out that the un- 
controlled smooth component in shell corrections, resulting from the incorrect treatment of 
continuum, can cause some renormalization of the parameters of the macroscopic energy 
formula, determined from fitting the nuclear masses. Unfortunately, the results discussed in 
Ref. ]S7| suggest that the renormalization procedure, based solely on bound states, produces 



very unstable shell corrections in cases when the upper limit of the averaging interval (i.e., 
the number of the highest single-particle level considered) approaches the actual number of 
particles. Actually, such a situation happens at the particle drip lines. 

The proper treatment of resonances is not an easy task, especially for deformed sys- 
tems. Therefore, other methods of dealing with continuum are more useful in practical 
applications. In 1973 Jennings proved |28j the equivalence between the Strutinsky approach 



and the semiclassical averaging based on the partition function method (Wigner-Kirkwood 
expansion). In the semiclassical approximation, the Fermi energy A sc is determined by 

N = £Z(^)> (2-19) 







and the smoothed energy of the system is 



E sc = NX SC - C£ (^P) , (2.20) 
where £ denotes the Laplace transform and Z is the semiclassical partition function [55] 

Z(f3) = ^J e"^ (l + h Wl + h 2 w 2 + ...) d 3 pd\ (2.21) 



6 



where H c is the classical one-body Hamiltonian of the system and the w's are defined in terms 
of the one-body potential. The partition function method was shown to be an excellent tool 
for computing the average single-particle energy [|30|-|33|j. In particular, Ref. |3(J contains the 



explicit numerical check of the equivalence between the Wigner-Kirkwood (WK) expansion 
and the full Strutinsky averaging with the correct treatment of resonances. (For more 
discussion, see Ref. [011, Table I.) 



III. WOODS-SAXON MODEL AND THE PARTICLE CONTINUUM 



The deformed shell model used in this study is assumed to be of the form of the average 
deformed Woods-Saxon (WS) potential, which contains a central part, a spin-orbit term, 
and a Coulomb potential for protons. All these terms depend explicitly on a set of external 
(axial) deformation parameters, (3\, defining the nuclear surface: 



R{9-(3) = C{l3)r A 1 l' i 



i + £/3 A y AO (0) 



where the coefficient C assures that the total volume enclosed by the surface (Eq. 
conserved. 



(3.1) 
1) ^ 



The deformed WS potential is assumed in the form [34 

V(r) = V [l + K I]f(r), 



(3.2a) 



™ = 1+ e* P W)A>r (3 - 2b) 

where the function £ (r) is equal to the (perpendicular) distance (taken with the minus sign 
inside the surface) between the point r and the nuclear surface represented by Eq. ( [3.1|) . 
The diffuseness parameter a is independent of Z and N. This choice of the WS potential 
guarantees that the diffuse region of the potential is constant independently of the nuclear 
deformation. 

The deformed spin-orbit (SO) potential is 

V S o(r) = -^(Vf so xp).s. (3.3) 
Here k is a (dimensionless) SO strength factor and 

/so(r) = — t (3.4) 

1 +exp [£so(r)/aso] 

is a WS form factor. (Since the radius of the SO potential, r D _so, is, in general, different 
from the radius of the central part, the distance function entering Eq. fl3.4|) is indicated by 
&<>•) 

The Coulomb potential for protons, Vc, is assumed to be that of the charge (Z — l)e 
distributed with the charge density pc( r ) = Pof(v), where f(r) is the WS form factor of 
Eq. ( p.2b| ). In our study we employed the set of WS parameters introduced in Ref. []35 
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These parameters have been widely used in nuclear structure calculations around the beta- 
stability line and for neutron-rich nuclei. So far, no attempt has been made to optimize 
these parameters for very neutron-rich nuclei. 
The eigenstates of the WS Hamiltonian 

H ws = T + V + V so + \ (1 + r 3 ) V c (3.5) 

are found by means of the expansion method in the axially deformed harmonic-oscillator 



basis. The oscillator frequencies, hu± and hu z , have been optimized according to Ref. |34 
In our calculations we used all basis states belonging to N<N osc (stretched) oscillator shells. 

The pattern of eigenstates of the WS Hamiltonian ( |3.5| ), as a function of the number 
of harmonic-oscillator quanta included in the basis, (10<iV osc <50) is displayed in Fig. 
The left portion shows the results representative of the spherical shape. [Note that every 
single-particle orbital is (2j + l)-fold degenerate.] The spherical symmetry is lifted in the 
right portion, representative of the deformed situation (/?2=0.4; here each orbital is only 
two-fold degenerate due to the time reversal symmetry). It is seen that the energies of 
bound states (e, < 0) converge very rapidly with the size of the basis; the convergence is 
achieved for iV osc ^14. The positive-energy quasibound states representing the discretized 
particle continuum vary dramatically with N osc . Asymptotically, as the basis becomes infi- 
nite, quasibound states approach zero energy (cf. Fig. 18 of Ref. |18j and related discussion). 
This leads to an increased single-particle level density at Cj > at large iV osc - values. The 
only states which are not strongly affected in the considered interval of N osc are the high-j 
orbitals which are fairly well localized inside a pocket in the centrifugal barrier created at 
positive energies by the mean-field potential. 

The situation in protons is slightly less dramatic than in neutrons. Because of the 
presence of the Coulomb barrier (~9 MeV in 120 Sn), the positive-energy single-particle states 
below ~5 MeV are fairly stable. (They represent narrow sub-barrier resonances, interesting 
in the context of proton emitters.) Judging from the results shown in Fig. ^|, it is impossible 
to find an "optimal" value of N osc , at which the positive-energy quasibound spectrum would 
become a fair representation of physical continuum (resonances). 

Since, at the particle drip lines, the radial asymptotics of wave functions are particularly 
important, it is instructive to relate N osc to the radial dimensions of the harmonic oscillator 
wave function. Its size depends both on the principal quantum number, N, and on the 
orbital quantum number, £. The classical major axis of the orbit (N, £) is given by |36 



\2 ^ 

max J ~ 



N + VN 2 - £ 2 } , (3.6) 



where L=^j^ is the oscillator length, i.e., for the £=0 oscillator states, the classical outer 

turning point is r mSLX ~L\ / 2N . A weak dependence of r max on N makes it very difficult, 
if not impossible, to describe the asymptotic radial behavior of the wave function using 
the expansion method in the harmonic oscillator basis. From this point of view, the direct 
integration of the Schrodinger equation in a spatial box is superior. For instance, for the 
nucleus 120 Sn discussed in Fig. 0, the oscillator length L appearing in Eq. (|3.6j ) equals to 
about 2fm, hence r max ~10fm for iV=12 and r max ~20fm for iV=50. (For a pedagogical 



example, see Fig. 19 of Ref. [18].) By the same token, for a high-£ orbit with £=N the major 



^7fm and 14 fm for iV=12 and N=50, respectively. 
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IV. WIGNER-KIRKWOOD SEMICLASSICAL EXPANSION 



The semiclassical approximation for a one-body Hamiltonian, including SO interaction, 
was developed in Ref. |30|| . Denoting by U(r) a one-body potential without the SO interac- 
tion, 



U(r) = V(r) + l(l + T 3 )V c (r) 
the particle number equation ( 2.19Q takes the form: 



(4.1) 



N 



l /2M\2 r» 



3tt 2 V h 



d 3 r{ (A sc - U)i + [| K 2 (V/so) 2 (Asc - U)'- 
- ^V 2 f/(A sc -f/)-5 



(4.2) 



The integral is cut off at the classical turning point (more precisely: the surface of turning 
points) defined by the relation U(r sc ) = X sc . The equation ( |4.2j ) for A sc is solved iteratively. 
The expression Q2.20 ) for the smoothed energy is slightly more complicated. In analogy to 
the notation introduced in Ref. PU|, we may write: 



Asr.iV 



E% + E\ + 25° J - [E s ? + E s ° 



where 
E\ 



157T 2 
1 



2MV 
IF) 



d 3 r(X sc - U) 



24tt 2 V h 



2M\2 r- 



1 



h 2 



d 3 r{X sc - U)^V 2 U, 
d r- 



(4.3) 

(4.4) 
(4.5) 



{-7V 4 £/(Vf/) 2 + 5(V 2 f/) ; 



5760tt 2 V2M; J (A sc -f/)t(W) 2 

+ io[w ■ v(v 2 u)]v 2 u - 5(v 2 u) 2 vu ■ v(yu) 2 /(yu) 2 

+ (V 2 U)V 2 (VU) 2 + VU ■ VV 2 (Vf/) 2 - V 2 (VU) 2 VU ■ V(VU) 2 /(VU) 2 } , (4.6) 



so 



E^ 
and 



67T 2 V h 



2M\2 rr, 



E 



so 
+1 



K 



h 2 



48tt 2 \2M) 
V/so • V(V 2 / SO 



d 3 r(X sc -U)HVf so ) 



d 3 r(X sc - f/)^Uv 2 (V/ SO ) 2 - (V 2 / SO ) 2 



(4.7) 



;V/so) 2 V 2 f/ 



2(A SC - U) 
|V/ S0 • V(V/so) 2 l + 2k 2 (V/so) 4 



2k 



;V/so) 2 V7so 



(4- 



Some technical details of the calculations of the above integrals, specific to our WS model, 
are discussed in Appendices A and B. 
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The WK expansion (( f4.2|) -(fO| ) ) is only applicable for one-body potentials whose first 
and higher derivatives exist. This means that the method cannot be directly applied to a 
wide class of single-particle models based on the average field obtained by means of folding 
over a sharp generating potential or a sharp density distribution. Examples here are the 
folded Yukawa potential |18[ or the Coulomb potential generated by the uniform density 
distribution; both are continuous in their values and first derivatives, but are discontinuous 
in their second and higher derivatives. 



V. RESULTS 

A. Quasibound states in the shell-correction method 

The dependence of the shell correction ( |2.10|) on the size of the oscillator basis and 
deformation is shown in Fig. |3] for a neutron-rich system ( 100 Zr) and a neutron drip-line 
nucleus ( 122 Zr). In the calculations all the quasibound states up to a cut-off energy of 
131.2 MeV/A 1 / 3 above the Fermi level were included. As expected, due to the effect of 
increased density of quasibound states with N osc , the shell correction is unstable. It decreases 
or increases linearly with N osc , depending on the position of the Fermi level. Indeed, for 100 Zr 
the value of -E^heii slowly increases with N osc , while the opposite effect is seen for 122 Zr. The 
reason for this different behavior is explained in Fig. ^ which shows the mean single-particle 
level density ( |2.9| ) as a function of single-particle energy e for 100 > 122 Zr. For relatively low 
values of N osc (e.g., A osc =12-14), the average level density increases monotonically and it 
depends weakly on p (results for p=6 and p=12 are very similar). On the other hand, for 
iV osc =50, there appears an artificial local minimum in g(e) located at -7<e<-5MeV which 
is caused by a rapid change in slope of g(e) due to a large number of quasibound states at 
e^O. This fluctuation is caused by a polynomial correction; as seen in Fig. |] the results for 
N osc =50 strongly depend on p. The smooth single-particle energy is obtained by integrating 
the product g(e)e up to the Fermi level A, cf. Eq. (|2.14|). Since <?(A)at osc= i4 < g(\) Noac=50 for 



122 Zr (opposite holds for 100 Zr), this explains the tendency seen in Fig. |3|. 

The instability of the shell correction with respect to iV osc means that the total mass of 
a drip-line nucleus cannot be estimated by means of the SAM. Moreover, since in different 
isotopes one obtains a different dependence on N osc , the two-neutron separation energies are 
also affected. In fact, as seen in Fig. [| the deformation energies can also be contaminated 
by the iV osc -dependence. The relative values of -E'sheii, computed at different deformations, 
vary with N osc up to iV osc =14 ( 100 Zr) and iV osc =18 ( 100 Zr), while the recommended value of 
A osc is 12-14 p|. 



The results presented in this section demonstrate that the uncontrolled error in -E s heii can 
seriously affect theoretical mass predictions for nuclei far from stability. If parameters of the 
macroscopic energy formula are determined by the global fit to experimental nuclear masses, 
some part of the unwanted effect of the continuum is absorbed by the isospin-dependent 
terms of the liquid-drop (or droplet) model mass formula. 



10 



B. Plateau condition in the shell-correction method 



In Fig. ^ the total spherical shell corrections for 100 Zr and 122 Zr (i.e., including proton 
and neutron contributions) are plotted as a function of 7 for different values of p. For the 
basis size we took the recommended value of iV osc =12. Neither for 100 Zr nor for 122 Zr is the 
plateau condition (|2.16|) fulfilled exactly. However, while in the former case the fluctuations 
of -E s heii with 7 and p are rather weak (for the "usual" values of p=6,8 and j/huj o =1.0-1.2 
the local plateau in E^n is at ~6.7MeV), in the latter case the variations in shell correction 
are more dramatic and no unambiguous values can be extracted. 

In the recent state-of-the-art large-scale mass calculations based on the folded Yukawa 
potential, Moller et al. [[| used r y=huj B s (B s is the ratio of the actual nuclear surface to 
the spherical surface, i.e., B s =l at spherical shape), p=8, and iV osc =i2. They noted that, 
in contrast to light nuclei which have very low density of single-particle levels, the plateau 
condition is fulfilled for heavy nuclei. Our analysis suggests that for nuclei far from stability 
the plateau condition can also be violated. 



C. Shell energies for the Wigner-Kirkwood smooth energy 

In order to estimate the continuum-related uncertainty in the shell energy in typical SAM 
calculations, we performed a semiclassical analysis using the WK expansion. Figure |6| shows 
neutron shell corrections E s ^ e n(SAM)=E s , p —E s , p computed by means of the SAM, Eq. ( |2.14| ) 



(A" osc =12, r y=hoj , p=8) as compared to E s h e \\(WK)=E sp — E sc where the smoothed WK 
energy ( [4.3|) is used. In both cases the same WS model was used. The calculations were 
carried out for the well-bound neutron systems 100 Zr and 166 > 186 0s (A n <-6MeV), and for 
neutron drip-line nuclei 120 ' 124 Zr and 206 > 226 Os (-4 <A n <-0.2 MeV). The average deviation 
between £ shcU (SAM) and £ shell (WK) ranges from ~0.4MeV in 100 Zr and 166 0s to ~4.6 WeV 
in 124 Zr and ~7.6 MeV in 226 0s. Generally, for the neutron-rich systems, the shell corrections 
obtained by means of the SAM are lower than the semiclassical ones and the difference 
increases with N. Therefore, in the standard shell-correction calculations the nuclei around 
the neutron drip line are overbound. Another interesting result presented in Fig. ^| is that 
the difference between £J s h e ii(SAM) and E'sheiiCWK) depends rather weakly on deformation. 
This suggests that one can approximate E s ^ e i\ by 

£ she ii(/3) = £ sh eii(/3; SAM) + £ shcll (/3 = 0; WK) - £ shcll (/3 = 0; SAM), (5.1) 

where the spherical semiclassical shell correction, £7 s h e ii(/3 = 0;WK), involves calculating 
one-dimensional radial integrals only. The renormalization ( |5.1| ) can be particularly useful 
in the large-scale global mass calculations, such as those of Ref. ||. However, one should 
keep in mind that this renormalization will not cure the problems related to the plateau 
condition and the dependence on N osc discussed above. 

It is worth noting that in the Extended-Thomas- Fermi-with-Strutinsky-integral approach 
|7|j37|-p9| the smooth energy entering their shell correction is defined by the semiclassical 



method ||38|| . Consequently, their shell corrections are free from the problems of the contin- 
uum. 
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D. Location of particle drip-lines in the shell-correction method 



In the macroscopic-microscopic method, the second contribution to the shell energy 
comes from pairing. Pairing correlations play a very special role in drip-line nuclei. This is 
seen from approximate HFB relations between the Fermi level, A, pairing gap, A, and the 
particle separation energies. For instance, the neutron separation energies S n and S^n are 



given by fl8|,f40 



S n ~ — A n — A, (5-2) 
S 2n « -2A n , (5.3) 

which leads to the following conditions for the one-neutron and two-neutron drip-lines: 

S n = A n + A = 0, (5.4) 
S 2n = A n = 0. (5.5) 



In particular, the condition ( |5.4j ) nicely illustrates the crucial role of pairing interaction for 
determining one-neutron drip line; it shows the equal importance of the single-particle field 
characterized by A (determined by the particle-hole component of the effective interaction) 
and the pairing field, A (determined by the particle-particle part of the effective interaction). 
In fact, just around the particle drip lines, particle-hole and particle-particle channels are 
very strongly coupled, and the standard procedure based on a "two-step" treatment of 
the pairing Hamiltonian (i.e., computing pairing correlations after determining the single- 
particle basis) seems inappropriate. 

In part II of our study fl(| , the pairing properties of drip-line nuclei are discussed using 
the self-consistent theory. Here, we only concentrate on the standard BCS or Lipkin-Nogami 
(LN) treatment of pairing correlations, which is used in the currently available large-scale 
mass calculations [|7|||. 

The macroscopic part of the total energy, -E macr , already contains the average pairing 
energy, which accounts for the main part of the even-odd mass difference. Therefore, it is 
the fluctuating part of the pairing energy 

^shell = -^pair ~ -Epair (5.6) 

that enters the expression for the total shell energy. The pairing correlation energy, impair, is 
usually computed using a monopole pairing interaction with a constant (state-independent) 
strength G. 

In this study, the BCS and LN equations were solved considering the N p (=Z or N) 
lowest single-particle orbitals. The average pairing strength, G avg , and the average pairing 
energy, -E pa i r , were calculated according to the average gap method [15,18|,4~i'l in the version 



of Ref. [I2"|. We have also performed calculations using the fitted pairing strength Gat- 

G at • A = G + Gi /, (5.7) 

where the constants G and G\ are obtained by a fit to odd-even mass differences in a given 
region of nuclei. 



12 



There is an obvious advantage of using G avg rather than G&t when going far from stability. 
Since G avg is inversely proportional to the average single-particle level density at the Fermi 
level, G avg oc g(A) -1 , this attenuates the influence of quasibound states on pairing properties. 
Figure [7] shows the neutron pairing gaps for the tin isotopes computed within the BCS 
approximation using G av g or G&f For neutron numbers 50<iV<84 both prescriptions for G 
yield fairly similar results. However, for higher neutron numbers the pairing gap obtained 
with G&t is much larger that that computed with G avg , and the result strongly depends on 
the size of the single-particle basis. In particular, for iV osc =40 the density of quasibound 
states is so large that A n ^0 for iV=82. Since the pairing correction (|5.6|) behaves as A 2 , the 



effect on E^^' is even more dramatic. Therefore, the contribution from quasibound states 
to the pairing correction is a source of another uncontrolled correction to the total energy 
in the macroscopic-microscopic method. 

The total energy of the shell-correction method, Eq. ( |2.1| ), is divided into three parts: 

E = -E macr + Ej^ + Ej£{* , (5.8) 



where E ma _ cr is the macroscopic energy of the Yukawa-plus-exponential model of Ref. |43| , |44 
E^cii is the shell correction (|2.10|) , and E^^ is the pairing correction computed using the 



BCS procedure. (Results of the LN calculations are fairly similar and will not be discussed 
here.) 

The predicted one- and two-neutron separation energies for the neutron-rich tin and lead 
isotopes are shown in Figs. [8] and || respectively. The two-neutron separation energies (top) 
were computed directly from binding energies of even-even nuclei. Since no self-consistent 
blocking was performed for odd nuclei, the one-neutron separation energies (bottom) were 
approximated by means of Eq. ( |5.2| ). As seen in Figs. |8| and |S], the influence of the continuum 
on the calculated positions of driplines is quite significant. The two-neutron separation 
energies calculated in the WK method are systematically lower than those obtained in the 
SAM and the difference approaches ~1.6 MeV at iV=106. This difference results in a shift 
in the position of the two-neutron drip line. According to the semiclassical approach, the 
nucleus 156 Sni 6 appears at the two- neutron drip line, while in the SAM calculations there 
are several (~5) more stable even-even tin isotopes expected. A similar trend is seen in the 
position of the one-neutron drip line which in the SAM calculations is overestimated by 2-4 
mass units. 



E. Fermi-level self-consistency condition 

For the lead isotopes the Fermi level in the WK method becomes positive at A r ~175, 
and consequently no solution to the particle number equation (|4.2j) can be found. At the 
same time, the value of S2 n for iV~174 is still positive (Fig. and equals to about 1.7 MeV. 
This constitutes a contradiction because, according to Eq. (|5.3| ), the values of Szn and A n 
should vanish simultaneously at the two-neutron drip line. 

In the self-consistent theory (e.g., HF+BCS or HFB), the Fermi energy is equal to the 
derivative of the total energy (ground state energy of the even-even system) with respect 
to the particle number. However, in the shell-correction method this relation is violated 
due to the particle-number inconsistency inherent to the macroscopic-microscopic model. 
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Differentiating both sides of Eq. (|5.8|) with respect to iV (assuming the parameters of the 
one-body potential to be fixed), one obtains five different Fermi energies, namely, 

Atot — A macr — A + A s . p . + SX, (5-9) 

where A tot is related to the neutron separation energy, A macr represents the macroscopic 
Fermi energy, A (or A sc in the WK method) is the smoothed neutron Fermi energy, A sp . 
is the neutron Fermi energy of the single-particle model obtained from the BCS equations, 
and 5X contains contributions from the smoothed pairing energy and the proton shell cor- 
rection. (The quantity 5X is small and consequently is ignored in the following discussion to 
simplify the presentation.) It should be noted here that the finite-difference approximation 
to derivatives of the single-particle energy with respect to particle number (i.e., A, A sc , and 
A s . p .) cannot be applied. (The deviation between derivatives and finite differences can be as 
large as lOMeV!). The reason is that in the model based on the average one-body poten- 
tial, the variation in the particle number (say from iV to iV+2) leads to (i) the variation 
in the chemical potential, and (ii) the variation in the average potential itself. Indeed, the 
parameters of the average potential depend smoothly on Z and N, and this causes a quite 
sizeable change in single-particle level density with particle number. 

In self-consistent approaches based on the two-body Hamiltonian, the requirement A to t = 
A s p . is fulfilled automatically. In the shell-correction method this requirement can be referred 
to as the Fermi-level self- consistency condition: 

A macr (Z, N) = X(Z, N) or A macr (Z, N) = X SC (Z, N) (5.10) 

for the SAM or WK methods, respectively. 

The parameters of the microscopic model are usually adjusted to selected single-particle 
properties of nuclei close to the beta-stability line, and the parameters of the macroscopic 
model are found by a global fit to masses and fission barriers. Therefore, it is not surprising 
that when extrapolating far from stability, the particle-number dependences of A macr and A 
are different. In Figs. [H] and [II] the Fermi energies defined in Eq. (|5.9|) are shown relative to 



Amacr for the tin and lead isotopes, respectively. The degree of inconsistency in A is measured 
by the magnitude of deviation AA = A — A macr . It is seen that |AA| is rather small around 
iV=82 and iV=126 for the tin and lead isotopes, respectively, and it reaches the value of 
about AA=1 MeV at drip lines. 

The condition (|5.10|) defines a certain coupling between the macroscopic and microscopic 
parameters in the shell-correction method. Since AA varies very smoothly as a function of 
iV (see Figs. [I(] and [II]), the condition (|5.10p it equivalent to 



Amacr -^min) X SC (Z, -/V mm ), A macr (^, iV max ) 0, (5.11) 

where, for a given atomic number Z, N m[n and iV max are neutron numbers corresponding to 
proton- and neutron-drip-line isotopes, respectively. (iV max can be computed from Eq. ( [4.2| ) 
after putting A sc = 0.) 

The consistency between the microscopic and macroscopic part of the energy in the shell- 
correction method has been discussed by several authors. In particular, Myers noticed HSJ . 
on the basis of the droplet model, that the parameters of the single-particle model (such as 
radius, potential depth, diffuseness) should be related to the parameters of the macroscopic 
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model. In the global mass calculations by Moller et al. the parameters defining the 
single-particle Hamiltonian do depend on the droplet model constants. In particular, as in 
Eq. (|5.16|) , the depth of the folded- Yukawa potential depends linearly on the average bulk 
nuclear asymmetry of the droplet model, 5. 

To illustrate condition (|5.11[) , we performed calculations based on the WS average po- 
tential with parameters of ref. [ 18| , adjusted according to Myers Namely, we used: 



In 5 r „ n 
«so = 0.9— fm 
In 9 



for the diffuseness, and 



with 



r -so 



Ro 



(a/R ) 2 I A 



1 1/3 



(5.12) 



(5.13) 



R = R P + 0.82 - 0.56/^p 
for the radius. The radius R p and the depth of the potential 



R p = 1.16A 1/3 (1 + e) [fm], 
V n , p = (52.5 =F 48.75) [MeV] 



depend on the relative neutron excess coefficient / through 

- 1(1 - 0.0056A 1 / 3 ) + 0.0028A 1 / 3 (1 + I 2 



1 + 3.15/AV3 
+ 0.335 2 + 0.00062A 2/3 (1 - I) 2 . 



0.147 



(5.14) 



(5.15) 
(5.16) 



(5.17) 



Figure shows AA=A SC — A macr for the lead isotopes. Different curves were obtained by 
multiplying by a factor b (=1.00, 1.05, 1.10, 1.15, 1.20) the neutron excess coefficient / enter- 
ing through 5 expression ( |5.16|) . (Direct renormalization of 5 is much less convenient since it 
leads to the modification of the isoscalar part of the potential depth.) As can be seen from 
Fig. [12|, the standard value, 6=1, gives the values of A sc too small, in particular when the 



neutron number approaches iV max ; the resulting value, iV max =194, is much larger than that 



of the macroscopic model, iV=182 (as indicated in Fig. [12]). However, by increasing the value 
of b by 10% conditions (|5.11| ) can be met for the considered lead isotopes. It remains to be 
investigated whether the Fermi-level self-consistency condition can be consistently fulfilled 
for other isotopes. 



VI. CONCLUSIONS 

The main objective of this study was to investigate the influence of the particle continuum 
on the shell-correction energies of the macroscopic-microscopic method. The shell corrections 
obtained by means of the standard averaging method were compared with those calculated 
with the semiclassical Wigner-Kirkwood expansion technique. The systematic error in -E^heii, 
due to the particle continuum, can be as large as several MeV at the neutron drip line. 
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According to our calculations, this error depends weakly on deformation. This suggests a 
possibility of renormalizing -E^heii only at the spherical shape. 

As demonstrated in our study, the use of quasibound states as physical resonances can 
lead to serious deviations when extrapolating off beta stability. In particular, the associated 
theoretical error in predicted particle drip lines can be as large as AA=10. We also em- 
phasized the role of the self-consistency condition between the microscopic and macroscopic 
Fermi energies. If this condition is violated, the relation between the Fermi energy and the 
separation energy is lost. 

We hope that our work will be helpful for future global calculations of nuclear masses in 
the framework of the one-body (macroscopic-microscopic) description. As will be discussed 
in the forthcoming study based on self-consistent two-body procedures [fToH , there are also 
other uncertainties related, e.g., to the choice of the effective interaction. A deep under- 
standing of the coupling between single-particle and pairing channels, and between discrete 
states and particle continuum, is a key to the physics of drip-line systems and a serious 
challenge for future work. 
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APPENDIX A: SEMICLASSICAL APPROXIMATION TO THE WOODS-SAXON 

MODEL 



The semiclassical equations (|4.2|) and ( f4.3|) are defined through the high-order deriva- 
tives of average potentials entering the WS Hamiltonian. In order to calculate VV^ and 
V 2 V, we took advantage of certain geometric relations specific to the definition of the WS 
potential (|3.2|). Figure |X3] shows the surface geometry typical for an axial system; here the 



nuclear surface, r=R(8), is given by Eq. ( pA\) . Given the radius-point r=r OB (see Fig. [13 
the function £(r) represents the distance ±|r^s| [taken with the minus (plus) sign inside 
(outside) the surface R(9)]. Denoting by n [n=V£(r)] and t the normal and tangent unit 
vectors to the surface at point A, respectively, one can write 

r = R + n£(r), (Al) 

where R=vqa- By acting with the gradient operator on both sides of Eq. ( All) , one obtains 



(Vn)((r) = 2 - VR. (A2) 

In order to calculate Vi?, one computes the gradient in the Frenet frame defined by three 
unit vectors n, t, and b=t x n. The result is 
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VR 



Q 



+ 



V 



where the curvature radius, g, and rj=\rAc\ are given by 



(A3) 



Q 



R 2 + R 



>2 



V 



R 2 + R 



R 2 + 2R' 2 - RR"' '' 1-R'cotO/R 
and R'=dR/de and R"=d 2 R/d6 2 . Finally, by means of Eqs. Q-i 



(A4) 



), one obtains 



Vn = V 2 £ 



1 



+ 



1 



(A5) 



In the spherical case, £(r) = r — R, g = 77 = R, and Vn = 2/r, as expected. For given 
r, the distance £(r) and the vector i? are calculated numerically by solving the equation 
d\r - R(6)\/d6=0. By denoting 



V {n) (r) 



d n V 



(A6) 



one can compute the derivatives appearing in Eqs 

W (n) = V {n+1) n, 
(yV {n) ) 2 = (v (n+1) ) 2 , 

y2y(n) = y(n+2) + y(n+l)y2^ 



and 



For example: 



(A7) 



The derivatives higher than 2 involve quantities such as n • V(V 2 £) or V 4 £ which are com- 
puted numerically by means of the finite difference method. It is known that the standard 
parameterizations of the nuclear shape give rise to a singularity of higher derivatives at r=0, 
or to their large variations at 8=0° or 0=90° H46,47j. In order to achieve better accuracy 
of higher-order terms of the WK expansion, it has ben assumed that the WS potential is 
constant for the values of £,(r)<(, cr i t = —7a. We have checked that up to /?2~0.7 the results 
are stable with respect to variations of £crit- 

It is interesting to note that although CVV) 2 =(V^) 2 is a simple analytic function of 
£(r), the Laplacian V 2 V depends both on £(r) and 9, see Eqs. ( |A7[ ) and This will 

lead to slight modifications of some discussion in Refs. p8| , |49| , where V 2 V was assumed to 
depend solely on the distance 

Special care should be taken when integrating singularities at A sc = U in Eqs. ( f4.2|) and 
( |4.3D . A practical way of handling singularities at at r sc is to employ the identity: 



d 3 t 



<$>(r) 



u 



+ 2 



d d r 

nin 

dtt$(r sc (tt),tt) 



3 $(r,ft)-$(r sc (fi),Q) 



u 



+ 21^^(0). 



\/A S c - t/(r min (fi),fi) 



r sc (n) r^X sc - U(r,Q) 
r mln (fi) (dU /dr) 



2 - 



r (d 2 U /dr 2 
(dU/dr) 



(A8) 
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where $(r) is a function of r = (r, Q) and the integration is performed over the volume 
surrounded by the surfaces r = r min (fi) and r = r sc (f2). In the above equation, the singularity 
at r sc has been removed at the expense of generating singularity inside the classical region 
at (V£/) 2 =0. (Such a situation happens for the protons, due to different radial behaviors 
of WS and Coulomb potentials.) In practice, however, the identity ( |A8| ) is only used in the 
narrow region around r sc where (Vt/) 2 is never zero. 



APPENDIX B: COULOMB POTENTIAL OF THE FERMI DISTRIBUTION 



The Coulomb potential generated by the WS (Fermi) charge distribution ( |3.2b| ) is given 

by 



V c (r) = I d 3 r' P ° - , w . , 1 , . (Bl) 

v ' J l + exp[f(r')/a] \r - r'\ K ' 

In the axial case this integral can be reduced to two dimensions (with integrand involving 
complete elliptic integrals). Since 

VVc(r) = -47rp c (r), (B2) 



the contributions from the Coulomb potential to higher-order derivatives in Eqs. ( f4.2| )- ([4.3| ) 
are calculated easily using expressions derived in Appendix 0. 

For a spherical shape, the angular integrations in Eq. ( |B~1| ) can be performed explicitly 
to give 

V c (r) = — dr'r' 2 pc(r') + 4vr / dr'r'p c (r'). (B3) 

r JO Jr 

Finally, the radial integrals are computed by means of identity 

1 / 1 + E^=i(-l) n exp[-n(i? - r)/a] if r < R, 



l + exp[(r-R)/a] \ ~ E^=i(-l) n exp[-n(r - R)/a] if r > R. 
The resulting fast- converging series gives Vc with desired accuracy. 



(B4) 
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FIGURES 



FIG. 1. Schematic illustration of coupling between bound states and particle continuum in 
drip-line nuclei. The potential wells are represented by the average Woods-Saxon field (plus 
Coulomb potential for protons). 



FIG. 2. Single-particle energies of the WS Hamiltonian ( |3.5|) for neutrons (top) and protons 
(bottom) in 120 Sn as functions of iV osc (the number of harmonic-oscillator quanta included in the 
basis). The solid (dashed) lines correspond to tt = + (ir = — ) orbitals. The left portion displays 
the results for spherical shape (/?2=0); here every single-particle orbital is (2j + l)-fold degenerate. 
The results characteristic of deformed shapes (/32=0.4) are shown in the right portion; each orbital 
is 2-fold degenerate due to time-reversal symmetry. 



FIG. 3. Neutron shell correction for 100 Zn (top) and 122 Zr (bottom) for three quadrupole 
deformations, /?2=0j 0.3, 0.6, as a function of the number of harmonic-oscillator quanta in- 
cluded in the basis. For the SAM calculations we have used 7=1.2 7lw =49.2 MeV/A 1 / 3 and an 
eighth order curvature correction (p=8). All single-particle levels lying below a cut-off energy of 
£c«t = 3.2 ?kj =131.2 MeV/^4 1 / 3 above the Fermi level were included. 



FIG. 4. Neutron average level density g for 100 Zn (top) and 122 Zr (bottom) for /32=0.3 as a 
function of single-particle energy e. For the SAM calculations we have used the value of 7=1.2 %uj Q 
and e cu t=3.2hto . The average Fermi energies A are indicated by stars. 



FIG. 5. Total shell correction at /32=0 as a function of the smoothing range 7 (in units of Hco ) 
for 100 Zn (top) and 122 Zr (bottom). The three different curves correspond to p=6, 10, and 14. The 
number of basis states used in the diagonalization is iV osc =12. 



FIG. 6. Neutron shell corrections of the WS model, 2£ s heii(SAM) and £J s h c ii(WK), given by 
Eqs. ( |2.14 ) and ( f4.3| ), respectively, as a function of quadrupole deformation for 100 > 120 > 124 Zr (top) 
and 166 ' 186 ' 206 : 226 Os (bottom). The SAM smoothing was performed with iV osc =12, j=huj , and 
p=8. 



FIG. 7. Neutron BCS pairing gaps for the even tin isotopes as a function of N computed using 
the average pairing strength G avg p2| (solid line) or the fitted pairing strength G& t , (Eq. (5.7) 
for G =19.2MeV and G\=— 7.4MeV, dashed line). The single-particle spectrum (including quasi- 
bound states) was obtained by diagonalizing the WS Hamiltonian in A osc =12 (top) and iV osc =40 
(bottom) oscillator shells. 
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FIG. 8. One- and two-neutron separation energies for the even neutron-rich tin isotopes 
predicted in the shell-correction method with the WS average potential and the Yukawa-plus- 
-exponential macroscopic model. The results based on the semiclassical Wigner-Kirkwood method 
(solid line) are compared with those obtained using the standard averaging including quasibound 
states. The smoothing was performed with iV osc =12 (dotted line), or iV osc =40 (dash-dotted line), 
7=^w , and p=8. 

FIG. 9. Same as in Fig. || except for the neutron-rich lead isotopes. 



FIG. 10. Fermi energies of Eq. (|5.9| ) relative to A 

macr <^S cl function of N for the neutron-rich tin 
isotopes. The degree of the inconsistency in A is measured by the magnitude of deviation A — A maC r 
(SAM variant) or A sc — A macr (WK variant). The particle numbers at which At o t=0, A=0, A sc =0, 
A m acr=0, and A s . p .=0 are indicated by stars. 



FIG. 11. Same as in Fig. 11 except for the neutron-rich lead isotopes. 



FIG. 12. Semiclassical Fermi energy, A sc , relative to A macr for the lead isotopes. The parameters 
of the WS model were taken according to ref. [|ll|. The curves correspond to different values of 
parameter b multiplying the neutron excess coefficient I entering through 5 the formula ( p. 16 ) 



for the neutron potential depth. The particle numbers N=N mstK at which A sc =0 and A macr =0 
[Eq. (|5.1l| )l are indicated by stars. 

FIG. 13. Geometric relations between a vector-point r, the distance function and the 

nuclear surface, R{6). The unit vectors n, t, and b (normal, tangent, and binormal vector fields of 
R) represent the Frenet frame field on R. 
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